### Sample size for lot quality assurance sampling

n.for.lqas <- function(p0, q=0, N=10000, alpha=.05, exact=FALSE){
  maxi <- nrow(cbind(p0,q,N))
  if(length(p0)==1 & maxi >1) p0 <- rep(p0, maxi)
  if(length(q)==1 & maxi >1) q <- rep(q, maxi)
  if(length(N)==1 & maxi >1) N <- rep(N, maxi)
  n <- N
  alpha.i <- alpha
  if(length(alpha)==1 & maxi > 1) alpha.i <- rep(alpha, maxi)
  if (exact) {
    # Hypergeometric distribution 2-by-2 table : See `help("Hypergeometric")'
    # 
    #   q   n-x      n
    # k-q   m-(k-q)  m
    #   k   N-k      N
    # where N = population
    #       n = sample size
    #       q = positive among sample
    #       k = total positive in the population
    m <- rep(1, maxi)
    k <- rep(1, maxi)
    for(i in 1:maxi){
      for(j in N[i]:1){
        m[i] <- N[i]-j
        k[i] <- trunc(p0[i]*N[i])
        if (dhyper(q[i], j, m[i], k[i]) > alpha.i[i]) break	
        n[i] <- j
      }
    }
    method <- "Exact"
  }
  else {
    # For normal approximation calculation
    # Formula: d=n*p0-z*sqrt(n*p0*(1-p0)*(N-n)/(N-1))
    
    for(i in 1:maxi){
      for (j in N[i]:1){
        if ((j*p0[i]-(qnorm(p=1-alpha.i[i]))*sqrt(j*p0[i]*(1-p0[i])*(N[i]-j)/(N[i]-1))- q[i]) < .001) break
        n[i] <- j
      }
      method <- "Normal approximation"
    }
  }
  n <- round(n)
  if(length(alpha) > 1 ) {alpha1 <- alpha }else {alpha1 <- NULL}
  if(length(N) > 1 ) {N1 <- N }else {N1 <- NULL}
  table1 <- cbind(p0, q, N, n, alpha1)
  colnames(table1)[colnames(table1)=="alpha1"] <- "alpha"
  table1 <- as.data.frame(table1)
  
  returns <- list(p0=p0, q=q, N=N, alpha=alpha, method=method, n=n, table=table1)
  class(returns) <- c("n.for.lqas","list")
  returns
}

### Print n.for.lqas
print.n.for.lqas <- function(x, ...){
  if(nrow(x$table) < 6){
    cat("\n")
    cat("     Lot quality assurance sampling","\n","\n")
    cat(c("                             Method =", x$method, "\n"))
    cat(c("                    Population size =", x$N,"\n"))
    cat("  Maximum defective sample accepted =", x$q, "\n")
    cat("     Probability of defect accepted =", x$p0,"\n")
    cat("                              Alpha =", x$alpha, "\n")
    cat(c("               Sample size required =", x$n, "\n","\n"))
  }else{
    cat("\n")
    cat("  Lot quality assurance sampling","\n")
    cat(c("  Method =", x$method, "\n"))
    if(length(table(x$N))==1) cat(c("  Population size =", unique(x$N),"\n"))
    if(length(x$alpha)==1) cat("  Alpha =", x$alpha, "\n")
    cat("\n")
    print(x$table)
  }}
